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merits in the network [p| ]l7| , |l8H the equations are written 

as, 
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It is shown that an arbitrarily weak (frozen) heterogeneity can induce global synchronized oscilla- 
tions in excitable media close to threshold. The work is carried out on networks of coupled van der 
Pol-FitzHugh-Nagumo oscillators. The result is shown to be robust against the presence of internal 
dynamical noise. 



Emergent synchronized oscillations is a key subject in a 
variety of fields ranging from physics to biology or med- 
ical sciences. In the last few years several papers have 
been published concerned with the possibility of trig- 
gering global oscillatory behavior through heterogene- 
ity and/or internal (dynamical) noise in the excitable 
medium Although these works are a significant 

step towards the understanding of emergent oscillatory 
behavior, many points remain unclear. For instance, in 
H a model was proposed to investigate emergent oscil- 
lations in pancreatic (3 cells in which half of the el- 
ements was in the silent phase while the other half was 
continuously active. Although the approach gave clean 
globally synchronized oscillations, it is doubtful whether 
such a symmetric arrangement may have any biological 
meaning. The combined effects of diversity and 
internal dynamical noise |)||) 10 have also been investi- 
gated |^,|6|,[ll],[l2]] . While, a common conclusion seems to 
be the significant role of dynamical noise in triggering 
global oscillations, a key point such as the size depen- 
dence of the results was not investigated in detail. Here 
we start from Cartwright approach jj| and explore the 
possibility of synchronization as a function of the amount 
of diversity (fraction of diverse elements). We show that, 
as the system approaches threshold for oscillatory behav- 
ior, the number of diverse elements required to trigger 
global oscillations becomes arbitrarily small. This is par- 
ticularly appealing from a biological point of view as the 
possibility of having a small number of cells different from 
the rest is always there. We also show that these results 
are not significantly affected by internal dynamical noise. 

We base our analysis upon the van der Pol-FitzHugh- 
Nagumo equations Jl3|-|l6| that, as discussed in M, are an 
adequate mathematical description of the circuit model 
(which involves a capacitor, a nonlinear resistance accross 
the capacitor, an inductance and a resistance) commonly 
used to represent a physiological excitable medium. In- 
cluding first-nearest-neighbors coupling between the ele- 
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where N is the number of elements in the network. All 
constants and variables are dimensionless. Variables ip 
and <f> are proportional to the potential accross the non- 
linear resistance (cell membrane) and the current through 
the supply, respectively. The subindex i indicates an el- 
ement in the network. The constant v is proportional 
to the potential supplied, (3 to the (membrane) resis- 
tance and 7 to the square root of the quotient induc- 
tance/capacitance. The coupling between the elements 
in the network is accounted for by k (see Q for a thor- 
ough discussion). In the model we allow constant v to be 
different on each element i of the network and to fluctuate 
dynamically (j)(t) is a gaussian noise and v§ a constant). 

In order to quantify the emergence of oscillatory behav- 
ior we calculate the spatiotemporal average of variable ip, 
namely, 



\ 



]^)EE^(*)-<^(*)> a ] 



(2) 



where < ipj(t) > represents the temporal average of the 
j-th potential. The initial time ti is chosen so that the 
contribution of transients to the average is minimized 
(note that for v near threshold the transient can be very 
long, see below) while tf is taken to cover a sufficiently 
large number of periods of the system in its oscillatory 
phase (for the values of the parameters given below the 
internal period varies around 10). 
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Synchronization was in its turn evaluated by calculat- 
ing the following average, 
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where site " 1" was randomly chosen. We could have ex- 
tended the sum to all pairs of elements < ij > but this 
would have prohibitively increased computation time in 
large networks. Note that using a s to test synchroniza- 
tion is far more demanding than most tests used in pre- 
vious analyses. In Eqs. (2) and (3) the discrete sum in t 
clearly accounts for (numerical) time averages. 

In the following we take /3=0.5 and 7=2 and vary 
the remaining parameters. In particular we investigate, 
for a given is, how a Q and <j s vary with the fraction x 
of elements with —v (hereafter referred to as impurity 
elements or, simply, impurities) distributed randomly 
in the network. In the absence of both internal noise 
(z/o=0) and coupling (k=0), oscillatory behavior occurs 
for \v\ < v c = 0.60412 (for the values of j3 and 7 chosen 
here, see ||). Calculations were carried out on L x L 
clusters (L — 10 — 40) with periodic boundary conditions 
and an integration step At = 0.002. 

In Figure 1 we plot the spatial average of ip for 10 x 10 
networks with v = 0.62 on all elements but two that have 
v = —0.62, with and without coupling (in the former case 
k = 0.5). First we note that, as remarked above, the sta- 
tionary state in the uncoupled case is only reached after 
rather long times (t > 100). Instead, in the coupled case 
the transient is very short and the system soon shows a 
coherent oscillatory behavior. Importantly, the existence 
of emergent oscillations do not depend on the actual lo- 
cation of impurity elements. In Figure 2 we show the 
average ip over the network and over five realizations of 
quenched disorder (different spatial configuration of im- 
purities). As the period of oscillation is weakly dependent 
on the location of impurity elements, the resulting pat- 
tern is a typical sum of oscillators with slightly different 
periods. It is interesting to note that even in the case that 
the two impurities lie at neighboring sites, global persis- 
tent oscillations emerge upon coupling. The results of 
Figure 1 are truly remarkable as global oscillations are 
promoted by a very weak diversity (2% in this case). 
Some characteristics of this central result are discussed 
in detail hereafter. 

Figures 3 and 4 show the parameters that character- 
ize the emergence of oscillation and synchronization (a 
and <7 S ), versus the fraction x of impurities, for v= 1.0 
and 0.61, i.e., far and close to threshold, respectively 
The results correspond to networks of linear size L=20 
and 40 with coupling constant k=2 and 8 for L=20 and 
40, respectively. This choice was motivated by the scal- 
ing argument of Ref. || (see also |ll]] and the discussion 
below), according to which one obtains solutions with 



similar properties in two systems of linear size L and 
aL if the diffusive coupling constant of the latter is in- 
creased by a factor a 2 (apart from border effects). Aver- 
ages were taken over five realizations (some checks with 
up to twenty realizations led to similar results) and in the 
time range £=200-600. First we discuss the results with- 
out dynamical noise. The critical impurity fraction x c 
(value of x at which a steeply increases) for z/=0.61 and 
1 approximately lies at x c «0.006 and 0.2, respectively, 
the results being almost independent of size, particularly 
in the former case, although the sharpness of the transi- 
tion to the oscillatory phase increases with the size of the 
network, as can be noted in Figure 4. On the other hand, 
x c shows no dependence on the coupling constant 
indicated by the results of Figure 3 and 4 and other data 
not shown in the Figures (this is so once k is beyond a 
critical value, see ||). In fact, x c can be derived, within 
a more than reasonable approximation, from a simple 
mean field approach, according to which the onset will 
take place when < v >= (1 — 2x)v equals v c . This leads 
to x c = 0.5(1 — v c jv), which gives £ c =0.006 and 0.198 
for z/=0.61 and 1.0, quite compatible with the numerical 
results of Figures 3 and 4. Note that the parameter that 
characterize synchronization a s is significantly smaller for 
v = 0.61 ]l9]| . A plausible explanation for this behavior 
is that as for v = 1.0 the transition occurs at much larger 
impurity concentrations, clustering is more probable, in- 
creasing the difficulty of synchronizing the whole system. 

Dynamical noise does not qualitatively change the re- 
sults discussed above. Figure 3 shows results for v = 0.61 
and vq = 1 |20|. The most noticeable (quantitative) 



changes are: i) At x=0, a a is higher than in the ab- 
sence of dynamical noise, although it is still not suffi- 
ciently large so as to consider the system being in its 
oscillatory phase, ii) consequently, the transition is less 
sharp, and, iii) synchronization is decreased (larger a s ). 
These results are in apparent contradiction with several 
analyses which indicate that dynamical noise increases 
oscillation and synchronization fl|6|JTll. However, this 
may be well due to the non-optimal noise level which is 
required for coherent resonance and to the more severe 
measures of oscillation and synchronization that we have 
adopted. Morever, we note that in those studies nothing 
was said about whether the effect survives as the size of 
the system increases. Preliminar results indicate that in 
fact it does not, in line with the rather small increase 
of a s that dynamical noise promotes near x=0 and the 
decrease in synchronization. In any case the main con- 
clusion of this analysis is that dynamical noise does not 
modify the previous result, that is, the dramatic effect 
that a small number of impurities has in systems near 
threshold. 

A final point concerns the effect of the coupling con- 
stant. Results for networks of L=20 and 40, ^=0.61 and 
x=0.02, are shown in Figure 5. It is noted that a reaches 
its maximum (constant) value for a coupling constant 
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that is significantly lower for the smaller network. In 
fact this occurs at k «0.45 and 1.7 for L=20 and 40, 
respectively. This is in accordance with the expected be- 
havior (derived from the diffusive character of the cou- 
pling term) discussed above. The results for the synchro- 
nization parameter are similar: to reach the same small 
levels of a s (less than 0.5, say) the coupling constant in 
L=40 should be 4 times larger. Despite the usefulness of 
the scaling trick in numerical computations, one should 
recall that in realistic systems the coupling between el- 
ements is typically intensive, that is independent of the 
system's size, and it is determined by intrinsic properties. 
Generally our simulations show that, when the coupling 
constant is kept fixed, the emergent oscillations and the 
degree of synchronization are less and less pronounced 
as the number of consitucnts is increased. This worsen- 
ing occurs either when these effects are induced solely 
by noise and when they are triggered by diversity, as 
discussed here. As far as experimental results are con- 
cerned, the basic point for the synchronous behavior to 
be observed is the strength of the effective coupling con- 
stant with respect to the number of elements. Another 
important feature is the type of interaction. Indeed, it 
is possible that non-diffusive couplings may lead more 
efficient mechanisms of synchronization. 

Summarizing, here we have discussed the possibility 
of triggering global oscillations in close-to-threshold ex- 
citable media through an arbitrarily weak heterogeneity. 
The work was carried out by assuming the existence of 
two possible types of elements in the network, one silent 
and another continuously active. The results clearly in- 
dicate that when the system is near threshold, global 
(synchronized) oscillations emerge for a small number of 
diverse elements. Dynamical noise does not alter this 
conclusion, although it may introduce some significant 
changes such as a decrease in synchronization. 
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FIG. 1. Variable tp as a function of time in a 10 x 10 net- 
work of van der Pol-FitzHugh-Nagumo elements with /3 = 0.5, 
7 = 2 and v — 0.62 for all elements but two with v = —0.62 
located at random in the network. For those values of /3 and 
7 the threshold for oscillatory behavior in an isolated element 
occurs at v c = 0.60412. The results correspond to no coupling 
between the elements k = 0.0 (thick continuous line) and for 
k = 0.5 (thin line). In the latter case the average of ip over 
the whole network is shown. 
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FIG. 2. Variable tjj as a function of time in networks with 
the same parameters of Fig. 1 and k = 0.5. The result rep- 
resents an average over the whole network and over five real- 
izations (each corresponding to a random spatial distribution 
of the two elements with v = —0.62). 
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FIG. 3. Filled symbols: Parameter used to quantify the 
emergence of oscillatory behavior, as defined in Eq. (2), in 
a heterogeneous excitable media described by Eq. (1) with 
v = 0.61 versus the fraction of elements x in the network 
with v = —0.61 (averages were done over 5 realizations of 
the disordered network). The numerical results correspond 
to networks of size 20 x 20 and 40 x 40 (symbol size pro- 
portional to the linear size of the network). The rest of the 
parameters in the van der Pol-FitzHugh-Nagumo medium are: 
P = 0.5, 7 = 2, and values of k discussed in the text. Empty 
symbols: Same for the parameter used to quantify synchro- 
nization, as defined in Eq. (3). Circles: without dynamical 
noise. Squares: with dynamical noise (vo = 1). The lines are 
guides to the eye. 



FIG. 4. Same as Figure 3 for v = ±1.0. Only results 
without dynamical noise are shown. 




FIG. 5. Parameters used to quantify the emergence of os- 
cillatory behavior (a a ) and synchronization (a s ) versus cou- 
pling k. The results correspond to a fraction of impurities 
of 0.02, v = 0.61, /3=0.5, 7 = 2 and bidimensional networks 
of linear size 20 (u empty circles and a s stars) and 40 (a a 
thick continuous line and a s dotted line) . No dynamical noise 
was included in the calculation. Averages were done over 5 
realizations of the disordered network. 
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